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Abstract 



A toggle switch consists of two genes that mutually repress each other. This regulatory motif is active during 
cell differentiation and is thought to act as a memory device, being able to choose and maintain cell fate 
decisions. Commonly, this switch has been modeled in a deterministic framework where transcription and 
translation are lumped together. In this description, bistability occurs for transcription factor cooperativity, 
while autoactivation leads to a tristable system with an additional undecided state. 

In this contribution, we study the stability and dynamics of a two-stage gene expression switch within a 
probabilistic framework inspired by the properties of the Pu/Gata toggle switch in myeloid progenitor cells. 
We focus on low mRNA numbers, high protein abundance and monomeric transcription factor binding. 
Contrary to the expectation from a deterministic description, this switch shows complex multi-attractor dy- 
namics without autoactivation and cooperativity. Most importantly, the four attractors of the system, which 
only emerge in a probabilistic two-stage description, can be identified with committed and primed states 
in cell differentiation. We first study the dynamics of the system and infer the mechanisms that move the 
system between attractors using both the quasi-potential and the probability flux of the system. Second, we 
show that the residence times of the system in one of the committed attractors are geometrically distributed. 
We derive an analytical expression for the parameter of the geometric distribution, therefore completely 
describing the statistics of the switching process and elucidate the influence of the system parameters on 
the residence time. Most importantly we flnd that the mean residence time increases linearly with the mean 
protein level. This scaling also holds for a one-stage scenario and for auto-activation. Finally, we study 
the implications of this distribution for the stability of a switch and discuss the influence of the stability 
on a specific cell differentiation mechanism. Our model explains lineage priming and proposes the need 
of either high protein numbers or long term modifications such as chromatin remodeling to achieve stable 
cell fate decisions. Notably we present a system with high protein abundance that nevertheless requires a 
probabilistic description to exhibit multistability, complex switching dynamics and lineage priming. 
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1 Introduction 



During differentiation, a cell and its progeny cascade through a number of lineage decisions from stem cells 
over progenitor cells to mature functional cells. Many decisions are assumed to be binary and realized by a 
toggle switch, a simple cellular memory device. This network module consists of two genes, inhibiting each 
other via mutual promoter binding. In each differentiating cell, one gene will eventually win this biomolccular 
battle, inhibiting the other gene and subsequently activating its lineage-determining downstream targets. In 
hematopoicsis, the generation of blood cells, a series of gene switches has been found to determine the 
differentiation path of hematopoietic stem cells and to direct the ratio of mature blood cells (1, 2). The most 
prominent example in this context is the mutual inhibition of Gata-1 and Pu.l, two transcription factors 
responsible for the development of erythroid and myeloid blood cells from common myeloid progenitors (3 5). 

Due to its importance in development, toggle switches are subject to both experimental and theoretical 
investigations (for a review, see (6)). Using a deterministic framework under the assumption of large molecule 
numbers. Cherry and Adler (7) discussed criteria for working switches. More specifically Rocder and Glauche 
(8), Huang ct al. (9) and Chickarmane et al. (10) used a simple deterministic model of the toggle switch 
based on ordinary differential equations in order to describe the Pu.l-Gata-1 switch in hematopoicsis. A 
comprehensive overview and comparison of the different deterministic toggle switch models is provided by 
Dufl^et al. (11). 

All these studies focus on the steady states of the switch and the parameter dependent bifurcations in a 
deterministic framework. However, protein variations of a differentiating cell influence the dynamics of the 
decision making process and lead to stochastic transitions between the two steady states. This randomness 
is induced by gene expression noise, which has been shown to be ubiquitous in biological systems due to low 
molecule numbers (12). Thus, the probabilistic frameworks developed in order to account for gene expression 
noise (see (13) for a review) have to be applied in order to understand fundamental aspects of toggle switch 
properties. 

Probabilistic models of the toggle switch account for low copy numbers and intrinsic fluctuations. In 
(14), the dynamics of an exclusive switch, where two genes share the same promoter, is discussed within a 
probabilistic framework. A comparison of simple switch circuitries is given in Warren and Ten Wolde (15). 
Contrary to deterministic models, transitions between the two macroscopic regimes where one of the two 
genes dominates are possible due to the inherently noisy gene transcription (16, 17), even without cooperative 
binding of transcription factors (18). More recent contributions focused on analytic descriptions (19, 20), the 
switching time between macroscopic regimes for different regulatory realizations (16, 21, 22) or parameter 
regimes (17), boundaries for the switching time (23), or delay effects (24). Notably, all of these approaches 
are based on a one-stage model of gene expression, where DNA is directly processed into hmctional proteins. 
However, it has been shown that the characteristics of protein noise strongly depend on the underlying 
expression model (25, 26). 

In this contribution, we abstract the regulatory details of the prominent myeloid Pu.l-Gata-1 mutual 
inhibition. Contrary to common belief, which advocates the lumping of the two stages of expression, we 
show that the inclusion of both mRNA and protein leads to an interesting change in system dynamics. The 
probabilistic two-stage description exhibits complex multi-attractor dynamics without autoactivation and 
cooperativity. Remarkably, a recent study reported low numbers of mRNAs in single murine blood cells: 
Warren et al. (27) found around 10 transcripts of the Pu.l gene per cell in common myeloid progenitors. 
Based on these findings we study a probabilistic description of a toggle switch with low mRNA numbers, high 
protein abundance and in accordance with the known role of Pu.l, monomeric transcription factor binding. 
We deliberately choose the simplest toggle switch model and neglect autoactivation due to our ignorance of 
the logic of activation and inhibition at the promoter. However our results can easily be extended and are 
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discussed for the case of dimeric regulation and exclusive autoactivation. 



2 Results 

2.1 A toggle switch based on a two-stage model of gene expression 

We describe the mutual inhibition of two genes, further on called A and B, using a two-stage model 
of gene expression (25, 26) with mutual inhibition being realized as DNA-protein binding (see Fig. 1). 
This kind of switch has been implemented in vivo by Gardner et al. (28). The model can be repre- 
sented as a set of biochemical reactions for A and B, respectively, and a set of reaction rates a, /3, etc.: 



DNAa - 

mRNAA 

mRNAA 



DNAa + hiRNAa 



mRNAA + ProteinA 



<5a. 



ProteiuA 

ProteinA + DNAb 
j^js^^bound JA^ ProteinA + DNAb 



DNAb - 
mRNAB 



DNAb + uiRNAb 



7b, 



/3b 



mRNAe ^ mRNAs + Proteine 
ProteinB 

ProteinB + DNAa 



Qj^^bound 



Qj^^bound ProteinB + DNAa 



(1) 
(2) 

(3) 
(4) 

(5) 

(6) 



Reactions 1 and 2 correspond to mRNA transcription from an unbound promoter and mRNA degradation, 
respectively. Reactions 3 and 4 resemble protein translation and degradation. The last two reactions 5 and 
6 describe the binding and unbinding of a protein to the antagonistic gene and thereby the transition from 
an active to an inactive promoter and vice versa. Bound DNA lacks the ability to be transcribed. We 
would like to emphasize that here and t~ are rates rather than times. Note that we assume monomeric 
transcription factor binding as the simplest of regulatory interaction (which has recently been shown to be 
able to induce bimodal gene expression (29)). Our system's topology is symmetric with regard to the two 
genes, and so are the two columns of reactions 1-6 upon the exchange of gene labels A and B. 

This model of gene expression is a highly simplified abstraction of the complex processes in the cell. 
Condensing transcription into a single biochemical reaction does not account for the various steps required to 
transcribe a gene, e.g. the assembly of the transcription initiation complex, unwinding of DNA or transition of 
the polymerase to elongation phase. Postprocessing and transport mechanisms are also neglected. However, 
simplified models of gene expression have successfully been applied to experimental data, supporting the 
validity of these simplifications (9, 30, 31). 

Most commonly one will study the properties of the system in a deterministic framework using ordinary 
differential equations (ODEs) that describe the time-evolution of species concentrations (7-10). The ODEs 
can directly be inferred from reactions 1-6 assuming mass action kinetics: 



2 



= Tg (1 - C^a) - T^dAUB 
^TOa = Q!A<^A - TA^a 

at 

— TiA = ^ArriA - SaUa 
at 

+ T^(l - de) - r^fiBriA 



_1 



de = (1 - de) - T^dBriA (7) 



dt 

d 

dt 



rriB = ctBdB - 7BmB 



riB = ^bJt^b - SbUb 

+ Tg (1 - C^a) - T^dATlB 



(8) 

(9) 



where d» is the abundance of unbound DNA*, m* is the abundance of mRNA* and n, is the abundance of 
Protein* for * e {A,B}. 

Bound DNA is expressed in terms of unbound DNA due to mass conservation. Solving Eq. 7, 8 and 9 at 
steady state by setting all time derivatives to zero yields two solutions, one being biologically irrelevant due 
to its negative species abundances. Given non-negative initial conditions the system will always converge 
towards the positive steady state solution (32), given by (see Supporting Material for details) 
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nA(«^)=nB(«^)=-^(l-r?) (11) 

dA^''^ = dB^''^ = (12) 
l + r] 



with T] = ^ '^"st- -^^^ parameters are positive and for simplicity assumed to be symmetric for players 

A and B (a = oa = ctB, ■ ■ ■)• 

We now assess the stability of the positive solution Eq. 10 - 12 using standard linear stability analysis. 
To reduce the complexity of our system for the stability analysis, wo apply a quasi steady state approxima- 
tion to the DNA binding/dissociation process (oJa = = 0), reducing the dimensionality of our system to 
four equations. We evaluate the corresponding Jacobian at the single positive solution and use the Hurwitz 
criterion to verify that all its eigenvalues have negative real part. We conclude that the system has one 
stable positive fixed point but we cannot analytically exclude the existence of limit cycles. However, inspec- 
tion of the system's phase portrait (see Fig. S4 in the Supporting Material) indicates that no limit cycles 
exist. Summarizing, we showed that the deterministic model has only one steady state solution and is thus 
monostablc. 

However, since the deterministic approach is only valid in the limit of large numbers, small molecule 
numbers of DNA, mRNA, and possibly proteins advocate a discrete probabilistic description of the toggle 
switch. We define the state of the system at time f as a vector x{t), where Xi{t) G No is the abundance of 
species i at time t. Note that the state space is discrete as opposed to the deterministic model. To emphasize 
this difference we use the uppercase notation Da, Db, Ma, Mb, Na, Nb for the number of molecules of the 
respective species. We can describe how the probability V{x, t) of being in a certain state x changes over 
time by using the master equation of the system (33) 

^(a;,i) = X] b^x^o'T^ix' ,t) - Wa;'x'P{x,t)] . 
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The first term considers transitions from states x' with rate w^^' to state x, while the second term accounts 
for transitions from x to all other possible states x' with transition rates Wrj-irj.. The transition rates, also 
called propensities (34) are determined by the reaction rates and the number of reagents of the corresponding 
reactions (see Supporting Material 2 for an explicit form of the master equation of the system). 

Even though the master equation describes the dynamics of the system more accurately than ODEs (most 
obviously for low particle numbers) it is still an approximation of cellular dynamics as it assumes spatial 
homogeneity inside a cell and does not account for time delays. Still, the protein distribution predicted by 
the master equation of a two-stage expression model was indeed observed experimentally (35), supporting 
the stochastic two-stage model. 

Since the master equation for the switch is analytically solvable only for a number of approximations 
(sec e.g. (36)) and not integrablc for large molecule abundances, we simulate the system trajectories using 
Gillespie's algorithm (37). Each trajectory follows the master equation, and the set of infinite trajectories 
constitutes the distribution that solves the master equation. To obtain appropriate parameters values for 
stochastic simulation, we delineate upper bounds for synthesis parameters from biophysical arguments and 
adapt degradation parameters to fit desired molecular levels. Table SI in the Supporting Materials lists 
all used parameter values. Additionally, the analysis below has been conducted for a second, differently 
motivated, parameter set (see Fig. S3 in Supporting Material) and yields qualitatively identical results. 

2.2 Dynamics and quasi-potential 

In this section we discuss the main features of the switch dynamics. Contrary to the deterministic model, 
time courses of the stochastic toggle switch model show multistable behavior (Fig. 2A). Given the parameters 
in Table SI our toggle switch can adopt different attractors: The two attractors where one player dominates 
the other (called S'a and S'b depending on which player dominates) are clearly visible in Fig. 2A. A careful 
inspection of the timecourse and the probability distribution in Fig. 2A shows that there also exist two 
intermediate attractors where protein numbers are similar (A^a — Ab ~ 0). These attractors are called S\ 
and 5*3 from now on. In the timecourses of the system (Fig. 2A) one observes that the system frequently 
switches between the dominating and the intermediate attractors. 

To get a deeper understanding of the complex dynamics of the system the notion of a quasi-potential can 
be used. The quasi-potential U of the system is calculated through the relation U{x) ~ — log T' (x), where 
■p{<^s)(x) is the steady state distribution of the system. The number of dimensions of the state space where 
the quasi-potential is defined equals the number of species in the system. Here the probability 7"^*'')(.t) 
of a state x in steady state is estimated from 15000 stochastic simulation runs obtained by the Stochkit 
software toolkit (38). In Fig. 2B the projection of the quasi-potential on the A'a — A^b, Ma — -^b plane is 
shown. The four attractors S'a, "S'b, S\ and Sg can be seen clearly in the quasi-potential of the system. The 
two attractors S'a and Sb appear as basins at the lower left and upper right corner of Fig. 2 whereas the 
intermediate attractors S\ and Sg are located at the center and are not well separated. The dominating 
attractors can easily be distinguished from the intermediate attractors via parameter dependent thresholds 
Xa,Xb in the protein dimension (see Supporting Material 4). 

Importantly, one has to keep in mind that the system considered is out of equilibrium and that the 
dynamics of a non-equilibrium system are not entirely determined by the gradient of the quasi-potential but 
by an additional curl flux stemming from the non-integrability of the system (39). As a consequence barrier 
heights in the quasi-potential do not necessarily correlate with the probability of crossing the barrier. 

To understand the dynamics of the switch in more detail we therefore consider for each state x in the 
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state space the outflux F{x) acting on the the system at this point (16). We calculate the outflux as: 

F{x)=V^^^\x)J2T^{y\x){y-x), 

V 

where the probability 'P{y\x) of state y succeeding state x and the probability V'^^^\x) are calculated 
from stochastic simulations. Note that the outflux is different from the concept of field lines used in phase 
portraits of ordinary differential equations. The outflux F(x) is plotted as small arrows in Fig. 2 (vectors 
are normalized and circles correspond the origin of the vectors) for all states x with 'P'^^^\x) > 2.5 • 10~^. 
This indicates where the system will move from the current state on average. Due to this outflux the system 
enters and leaves the attractors Sa and Sb through different paths. This phenomenon has been described in 
(40) and linked to the emergence of time directionality in non-equilibrium systems. In order to move from 
high {Sa or Sb) to low (So) protein numbers, at first the corresponding mRNA number has to drop. On the 
contrary, moving from low to high protein numbers requires the rise of mRNA numbers first. 

A different view on the system's dynamics is provided by the quasi-potential landscape and outflux in the 
^totai^ ^totai pjj^j^p ^pjg where iV*^°tai = (i _ _^ is the total number of ProteiuA in the system, 
bound to DNA (first term) or free (second term). Choosing iy*°*^' and iVg°*^' as projected dimensions shows 
four distinct basins in the quasi-potential landscape. Two basins correspond to the attractors 5a and 5b. 
These are characterized by high amounts of the dominating protein and zero proteins of the repressed species. 
The attractors 5^ and 5g are now clearly separated. In these two basins a single protein of one species is 
present and only a moderate protein number of the other species. In the following we show why these basins 
emerge and how the system moves between the attractors. 

We explain the dynamics of the system with a typical trajectory of the system: Let us start with 
the trajectory in the attractor 5a (lower right) where ProteiuA dominates Proteine. Due to stochastic 
fluctuations in the promoter status, eventually a burst of proteins of B will occur and inhibit the promoter 
of A, whose protein numbers will drop (Fig. 3, trajectory I). While the formerly dominating ProteiuA's are 
degraded, also the newly created Proteine quickly decreases in numbers and only one bound Proteine is 
saved from degradation. This drives the system towards the origin in the quasi-potential of Fig. 3. However, a 
single ProteiuB cannot completely suppress the promoter of DNAa, leading to a small but constant synthesis 
of ProteiuA. The system settles into an intermediate state (5a) defined by the presence of one Proteins and 
an intermediate amount of ProteiuA originating from the leaky inhibition of DNAa and bursting. In order 
to leave this basin the system has two options: Either the single Proteiue is degraded when it momentarily 
is not bound to the promoter. Consequently the levels of ProteiuA rise again and the system reaches 5a. 
The system is moved to the lower border of the quasi-potential where a strong outflux pushes it towards 5a 
(Fig. 3, trajectory II). Alternatively, a burst of ProteiuB displaces the system from 5a into regions where 
the vector field points strongly towards the diagonal A^a*"*''' = ATb*"*^' (Fig. 3, trajectory III). However this 
burst is typically not strong enough to move the system onto the diagonal and it will fall back into the basin 
5^. In order to enable a change from 5^ to S^ the system has to reach the diagonal. This is accomplished 
if, while the system is moving towards the diagonal after the burst, additional bursts of ProteiuB move it 
onto the diagonal (Fig. 3, trajectory IV). Once the system has hit the diagonal both protein levels will drop 
to very low numbers since non of the players has any significant advantage. Here by chance the system will 
move to any side of the diagonal and either towards 5a or 5^ (Fig. 3, trajectories V, VI). 

We find that leaving 5^ towards 5a (Fig. 3, trajectory II) is much more probable than hitting the diagonal 
from 5^ (Fig. 3, trajectory IV), which would provide the chance of switching. This is obvious from the 
mechanism described above: Even though the events triggering the two alternatives (degradation of ProteiuB 
and an initial burst of ProteiuB) have similar probabilities, the diagonal crossing requires additional events 
and is therefore much less probable. This cannot be deduced from the quasi-potential landscape alone: From 
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Fig. 3 it can visually be inferred that the barrier separating and is higher than the barrier separating 
5^ and S^. This wrongly suggests that moving between 5^ and occurs more frequently than moving 
between Sa and 5^. 

Comparing the system dynamics of our switch with other descriptions we find that (i) deterministic one- 
stage and two-stage models show no bistability while (ii) a probabilistic one-stage model exhibits tristability 
with only one intermediate attractor (see Fig. S2 in the Supporting Material and (29)). We speculate that 
translational bursting destabilizes the intermediate attractor of the one-stage model, where none of the two 
players can overwhelm the other. Bursting provides an easy mechanism to escape this deadlock situation: 
It gives the player whichever bursts first a huge advantage over the other, giving rise not only to one 
protein (as in the one-stage model) but several proteins. As a result, the two-stage system is always quickly 
pushed away from the diagonal and stabilizes in the attractors 5^ or S^. Thus, only the combination of 
a probabilistic description with a two-stage model of gene expression leads to the complex multi-attractor 
dynamics described above. 

2.3 Residence times 

Genetic toggle switches are thought to be involved in the differentiation process of cells. A common idea is 
that different cell fates correspond to the different attractors of the system (41). Therefore it is of interest 
how long the system will stay in one of these attractors. In this contribution, we focus on the time the 
system will stay in the attractors Sa or Sb- We assume that only in these two attractors the concentration 
of either player is sufficiently high to carry out a downstream biological function which resembles the switch's 
decision. 

In previous contributions, such quantities have been calculated or determined by stochastic simulation 
for simpler switch models and were called spontaneous switching time (23), switch lifetime (15), mean first- 
passage time (14), or switching time (22). Since the switch may fiip from a dominating to an intermediate 
attractor, we choose residence time as the appropriate term for the quantity calculated below. In the 
following, we derive an analytical approximation for the time the switch stays in a dominating attractor, Sa 
or Sb, called the residence time tg. A simulation study for S]^/S^ suggests qualitatively similar behavior 
(see Fig. S5 in the Supporting Material). 

Let us assume that the system is in attractor Sa- Hence, the promoter of DNAb is bound by ProteinA 
while the promoter of DNAa is unbound. We assume that the protein levels in this attractor can be 
described with the simple two-stage model (25), resulting in a mean ProteinA level of Na = (qa/3a)/(7a^a)- 
Consequently, the protein level of Proteins is A^b = as it is inhibited by the high levels of ProteinA. 
In order to leave Sa it is crucial that one Proteine is synthesized, which then can bind the promoter of 
DNAa and shut down the synthesis of ProteinA, ultimately driving the system out of Sa and into S"^. This 
trajectory (called trajectory I in Fig. 3) involves the following events: (i) unbinding of ProteinA from DNAb, 
(ii) synthesis of ProteiuB during the unbound phase, and (iii) binding of ProteiuB to the promoter of DNAa 
before ProteiuB is degraded. 

First we describe the unbinding of ProteinA from DNAb. While the system is in Sa, ProteinA dissociates 
various times, leaving the promoter of DNAb unbound. The average time the promoter remains unbound, 
tu, is equal to the average time until a binding reaction occurs, which is 

The time the promoter stays unbound is a random variable itself, but for simplicity we approximate it with 

its mean value. Note that t„ depends, somewhat counterintuitively, on t"*" and not on t~, with t^^Na being 
the propensity for a binding reaction. Again we emphasize that and r" are rates (rather than times). 
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To ultimately synthesize a Proteine, at least one mRNAB has to be transcribed during t„ and translated 
before degradation. The probability of k transcription reactions to happen during t„ is 

-Ppoisson(-fS^ = k) = ^"'^ jj""^ ■ exp(-Q!B • tu) , 

as the number of transcription reactions K during tu is Poisson- distributed with mean ae • tu- Thus, the 
probability of at least one transcription during the unbound phase is 



qs = l- P{K = 0) = 1 - exp 



aB 



The probability of translation during an average mRNA lifetime I/7B is accordingly = 1 — exp (— /3b/7b)- 
Finally the probability for a binding reaction during average protein lifetime 1/Sb is qb = 1 — exp {—t^/Sb) ■ 
However, not only one but several unbound phases may occur before ProteinB is successfully synthesized. 
The number L of unbound phases until and including successful synthesis follows a geometric distribution, 
P{L = Z) = (1 — qy~^q with parameter q = qg-q^- q^. The average number of unbound phases during a time 
interval Ai is - At. Thus, we can convert the random variable L into T = L/tJ^ via a linear transformation 
of a random variable, giving the actual time until successful synthesis of ProteiuB. Notably, the derivation of 
the distribution for residence times goes beyond previous mean-field approximations. Using the properties 
of the geometric distribution for the random variable T, we end up with the mean and the variance of the 
residence time: 

, 1 ^ 2 1 i^-qsqtqb .,„^ 
ts = — and cTj = — — • . (13) 

ta • qsqtqb {Tp^ )2 {qsqtqbY 

An important approximation for the residence time can be derived under the assumption of rapid trans- 
lation and slow mRNA degradation, /3 ^ 7, leading to w 1. This implies that it is quite certain that 
an mRNA will be translated at least once before degradation. In the regime of rapid transcription factor 
binding (r+ » 5, a) the probability for a binding reaction is close to one, « 1, while the probability for 
at least one transcription can be approximated with qg ~ q:b/(ta-^a)- Taken together, this leads to a linear 
dependence of the residence time on the protein number, 

ts « (t+Zta) • (A^A/ae) • (14) 

We want to compare our analytical approximation with the residence time derived from simulations. To 
that end, we have to infer the dominating attractors from the simulated time courses. Recall that we can 
identify the dominating attractors via thresholds xa,Xb at protein levels. The residence time of attractor 
5'a {Sb) is estimated as the consecutive time in a trajectory where Nx > XA {Nb > xb)- Wc compare the 
analytically derived geometric distribution for the residence times (see Eq. 13) with numerical results by 
simulating the switch with a given parameter set and estimating the residence times from 10000 stochastic 
simulations. Fig. 4A shows excellent agreement between the geometric distributed residence time and the 
simulations for a protein degradation rate of (5 = 8 • 10~^s~^. This legitimates the approximations and 
assumptions made above for the parameter regime of rapid transcription factor binding. From the analysis 
of the mean residence time for different protein half-lives, we find again a good agreement between the 
simulation and the approximation (see Fig. 4B). Moreover, the slope of the log-log curve of the simulation 
is 1 — confirming a linear dependence of the residence time from the mean protein level. 

With the result from Eq. 13 we can compare the mean residence time of different switch models. First 
we consider a gene expression model where transcription and translation are condensed into a single protein 
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synthesis reaction. In analogy to the two-stage model of gene expression (26), this can be called a one- 
stage model of gene expression. To achieve the same amount of proteins at similar degradation rates, the 
synthesis rate in the one-stage model needs to be larger compared to the transcription and translation rates 
in the two-stage model. The probability qt which accounts for translation during mRNA lifetime can be 
set to 1, since there is no mRNA stage and proteins arc produced immediately. The binding probability qh 
remains unchanged. However, because of the increased synthesis rate, the probability qg of synthesis during 
the unbound phase will be larger than in the two-stage model. Therefore, the mean residence time will be 
decreased in the one-stage model as compared to the two-stage model, leading to more frequent attractor 
changes. This finding is in accordance with the previously reported stabilizing effect of bursts in an exclusive 
switch (16). 

A second modification of the switch includes not only miitual inhibition but also aiitoactivation of both 
genes. If the promoter of the gene is unbound it will be transcribed with a small basal rate k. If the promoter 
is bound by its own protein product the gene will be transcribed with full rate a ^ k. Repressor bound 
promoters arc inactive. For simplicity wc assume that cither activators or repressors are bound but not both 
at the same time. Note that in this case also the deterministic ODE model is bistable (42). Considering 
the mean residence time in a two-stage switch with autoactivation, we find that the probability qg of mRNA 
synthesis during the unbound phase is smaller than in the ordinary two-stage model. Since no activator is 
present in this attractor, mRNA has to be transcribed with the small basal rate k, making the transcription 
more improbable. The probability qt for translation remains unchanged. However, the probability qb of 
protein binding to the antagonistic promoter is also decreased since this promoter is occupied by the abundant 
activator most of the time. Therefore, repressor binding to this promoter requires an additional dissociation 
reaction of the activator during repressor lifetime. As both and qb are decreased the mean residence time 
in switch models with autoactivation will be strongly increased compared to the ordinary two-stage model. 

Summarizing, we find that the residence time is (i) geometrically distributed, (ii) the mean of the distribu- 
tion grows linearly with the number of proteins for slow mRNA degradation, and (iii) both the intermediate 
step of mRNA production and the autoactivation of transcription factors increase the residence time. 

3 Discussion 

3.1 Lineage priming 

We now discuss the implications of our findings in the context of cell differentiation driven by the toggle 
switch. In previous studies (8, 9, 40), attractors where either one or the other player is dominating, thereby 
repressing the antagonist {Sa, Sb) corresponded to committed cells. We also find analogs for the intermediate 
states S*^ and 5*3. In these attractors the system has a strong preference towards one specific dominating 
attractor, but is not fully committed yet. A similar behavior is known as lineage priming in stem cell biology 
(43). Two different studies (44, 45) showed that a population of stem and progenitor cells, respectively, can 
be divided into subpopulations that mainly give rise to only one of two possible cell types. In our simple 
model this would correspond to stem cells that reside either in 5^ or S^. These stem cells can still give rise 
to both cell fates but have a strong tendency towards one of them. 

Remarkably only a two-stage probabilistic model of the toggle switch shows dynamics reminiscent of 
lineage priming. Although a progenitor state exists in one-stage models of the toggle switch, cells in this 
state will move to either the one or the other committed state with equal probability. 
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3.2 Residence time 



We find that the residence time in 5a and 5b, a key property of the system, is geometrically distributed. 
Previous contributions (22, 23, 46) focused only on the mean residence time and did not consider its imdcr- 
lying distribution. What does a geometric distribution for the residence time imply for the differentiation 
process dependent on the state of a genetic switch? To discuss this question, let us first reason on how a 
differentiation decision could be established with the toggle switch lined out in the previous sections. 

We discriminate two scenarios for the differentiation of a cell: In the first scenario, the state of the switch 
completely determines the cell fate. Starting in the progenitor attractors 5^ or 5^, after a certain amount 
of time, the switch will move to a committed attractor. We assume that the high numbers of proteins of the 
dominating player will trigger the differentiation program of the associated lineage and establish the mature 
cell type. However due to stochasticity, the switch will drop out of the committed attractors and the cell 
will not only lose the current lineage decision, but possibly even switch to the opposing cell fate. In order 
to establish stable lineages in this scenario, the cell has to assure that the residence time of the switch is 
much longer than all relevant biological processes of the cell, especially cell lifetime. This guarantees that 
the cell will keep a lineage decision once it has obtained one. Yet the geometric distribution of the residence 
time imposes difficulties in this scenario: Even if the mean residence time is high, short residence times 
will always be more probable than longer residence times. The toggle switch could either be stabilized with 
the aforementioned autoactivation of the players, or with very high protein numbers so that the g(X)nietric 
distribution flattens and transforms to an almost uniform distribution. Both means would assure that only 
a very small percentage of a population of cells forgets its lineage decision during lifetime. 

In the second scenario we assume that the cell gets locked in one fate by changing the shape of the 
underlying potential so that further transitions between attractors become less possible. Such a change of 
the potential could for example be facilitated by chromatin changes, as proposed by Akashi et al. (47). In the 
following, we assume that only if one state dominates the other for a long enough fixation time td, downstream 
genes necessary for the decision process are activated (e.g. leading to chromatin remodeling), and the cell 
differentiates. Such a time depending property could be implemented with low-pass filters (see (48) for 
an example in hematopoietic stem cells) and would allow for an integration of external signals (see (49) 
for the instructive power of hematopoietic cytokines). In this scenario, the residence time determines when 
differentiation will occur: The switch will constantly move into and out of the dominating attractors, until the 
residence time is finally long enough so that the dominating player can activate the downstream differentiation 
machinery. Ignoring the time the system spends in the intermediate attractors and just summing up the 
residence times in 5a and 5b until a long enough residence time for differentiation occurs, we find that this 
time follows a geometric- like distribution (see Supporting Material 5) . Under this differentiation mechanism, 
most cells will differentiate very fast and only few cells will need longer. Experiments that measure the time 
for single cells needed to go from the primed to the committed state (as an extension to the 2-day threshold 
reported by Heyworth et al. (50) for GM-CTC cells) in order to support or reject these hypotheses remain 
to be done. 

3.3 Compcirison to previous models 

Finally we discuss how our findings relate to previous studies on the toggle switch. We found that the 
mean of the residence time distribution scales linearly with the mimber of proteins in the system. The more 
proteins are present, the longer the average residence time in 5a or 5b. However, shorter residence times 
are still most probable due to the geometric distribution. This holds for the one-stage, the two-stage, and 
the auto-activating scenario. 

This linear scaling differs from the exponential (23) or near exponential (46) scaling described previously 
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in the one-stage scenario. In contrast to our model, the model of Warren and ten Wolde (46) considers 
dimerization of the transcription factors, motivated by the fact that cooperative binding is necessary to 
achieve bistability in a deterministic framework (10). We showed that, as soon as stochastic fluctuations are 
introduced, a system with multiple attractors is achieved that can act as a proper switch with additional 
states of low co-cxpression. Including dimerization as a prerequisite for inhibition in a one-stage model will 
strongly increase the stability of the attractors S'a/b (46). This is consistent with our findings: Instead of 
requiring translation of one protein of the suppressed species, we now require this rare event to happen twice 
during a short time frame, which is much less probable. However, the inclusion of dimerization will have less 
effect on the two-stage switch: Since proteins are typically synthesized in bursts (in our model the average 
burst size is /3/7 = 10) and dimerization is a fast process (46), as soon as one burst occurs almost certainly 
a dimer is formed and can inhibit the currently dominating player. Therefore the probability of leaving the 
attractors Sa/b is similar to a non-dimeric inhibition. 

Contrary to our results, Warren and ten Wolde (46) report that introduction of mRNAs reduces the 
stability of the switch. This discrepancy can be understood in the light of dimerization. In their one-stage 
model dimerization is a key ingredient of stability, which is lost when introducing translational bursts ("shot 
noise"). As we considered monomeric transcription factor binding, stability does not rely on dimerization. 
Therefore mRNAs increase the stability of the system, because they introduce additional conditions required 
for switching. 

Due to these differences in the model it is hard to resolve the discrepancy between our linear and the 
exponential scaling of residence time foimd by Bialek (23) and Warren and ten Wolde (46). However we 
want to emphasize that the theoretical results shown in (46) only consider protein numbers up to 30. In this 
region our simulation results show slight deviations from the analytical linear dependence (Fig. 4). At such 
low protein numbers the system does not only leave the dominating attractor according to the mechanism 
described in our results. It is also likely that just due to fiuctuations in the gene expression (not fluctuations 
in the promoter) the dominating attractor is left. This mechanism operates only at very small protein 
numbers and its probability rapidly decreases with rising protein numbers. Therefore our results do not 
contradict the findings of Warren and ten Wolde (46), but consider a different parameter regime with higher 
protein numbers. Interestingly, the noise-driven attractor changes are also described by Kashiwagi et al. 
(51) where the authors link this mechanism to the selection of a favorable, less noisy attractor in E. Coli 
populations. 

In another contribution, Morelli et al. (52) use the forward flux sampling algorithm to assess the stability 

of a one-stage genetic toggle switch with dimeric transcription factor binding. They find a similar mechanism 
of attractor flipping which is based on the synthesis of the suppressed species due to promoter fluctuations. 
Using the forward flux sampling, they obtain estimates of the switching rate (the inverse of the mean 
residence time) for different amounts of fluctuations in DNA-protein interaction and dimerization. Morelli 
et al. (52) modulate the size of fluctuations at the promoter by varying the ratio of binding rate and 
synthesis rate, the adiabaticity parameter w = t'^ /a (t~ is adjusted to keep t'^ /t~ constant). Small to 
leads to strong fluctuations, whereas large lo reduces fluctuations. They find that increasing lo decreases the 
average switching rate and therefore stabilizes the switch. This dependency vanishes for w > 5, where the 
average switching rate remains constant. The latter is in accordance with our results in Eq. 14, where the 
mean residence time depends only on the ratio of r"*" and t~, not on the absolute values and is therefore 
independent of w. The dependency of the average switching rate for w < 5 is not predicted by Eq. 13 and 
14. It is also not visible in the stochastic simulations, where mean residence times of systems with lj = 1 and 
a; = 20 coincide (Fig. 4). The results of Morelli et al. (52) were simulated for an average number of proteins 
N\ = Nb = 27. As mentioned above, in regions of very small protein numbers the system might leave the 
dominating attractor by a mechanism not captured by Eq. 13 and 14, probably causing the difference of the 
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results of Morelli et al. (52) and our results for small u). 
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Figure 1: Scheme of the two-stage switch. Species associated with gene A arc shown in white, species 
associated with B are shown in gray. Sohd arrows indicate synthesis and binding, jagged arrows indicate 
degradation. mRNAA is transcribed from DNAa with rate aA- It decays with rate 7a and is translated 
into ProteinA with rate /3a- ProteinA decays with rate Sa and can bind (unbind) DNAb with rate tJ^ (^a )• 
Protein-bound DNA leads to transcriptional arrest. The topology is symmetric with respect to the genes A 
and B, thus, the same reactions exist for B. 



15 




M -M„ 

A B 

Figure 2: Dynamics and quasi-potential of the switch showing the different attractors of the system. (A) 
The timecourse of NA{t) — N^it) clearly shows the dominating attractors, which can be separated in state 
space via the thresholds xa and xb- Either A dominates (attractor Sa), or B dominates (attractor Sb), 
or the system is temporarily locked by two bound promoters with only marginal protein expression of A or 
B (attractors 5a and S^). A histogram of Nji^{t) — N^it) is shown on the right. (B) The quasi-potential, 
defined as U{x) = —logV^'^^\x), includes the mRNA dimension of the system. It shows the four possible 
attractors as basins in a probability landscape. Sa and Sb are visible as basins at the lower left and upper 
right corners, whereas and 83 are located around the origin ( Na — Nb = Ma — Mb = 0) of the landscape. 
Additionally, the outflux F{x) acting on the system at the state x in state space arc indicated as lines (circles 
correspond the the origin of the vector). Note that the outflux is different from concept of deterministic field 
lines. These vectors show that there are different paths for entering and leaving the dominating attractors. 
Parameters for the simulation are given in Table SI. 
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Figure 3: Quasi-potential of the system projected onto the Nj^ and iVg"* dimensions. Note that both 
axis are on logarithmic scale and are shifted by 1 in order to include A^X"*""' = and 7V^°*^' = 0. Therefore 
the lowest row in the plot corresponds to the case iVg°'''' = 0. The quasi-potential U{x) — —logV''^^\x) 
is color coded where red areas reflect minima of the landscape. Visible are four minima corresponding to 
(lower right), Sb (upper left), S'^ (lower middle) and (middle left). The vectors of the outflux at 
each point in state space are drawn as lines (circles correspond the the origin of the vector) . Note that the 
outflux is different from concept of deterministic field lines. In contrast to Fig. 2 the vectors are normalized 
and therefore show only the direction, not the magnitude of the field. Bold arrows reflect typical trajectories 
(I- VI) of the system. For a discussion, see the main text. 
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ts in hours Protein level 

Figure 4: Residence time tg in the two-stage toggle switch. (A) The distribution for tg obtained by stochastic 
simulation is in good agreement with the geometric distribution derived from our mean-field approximation. 
The mean of the distribution is indicated by a dashed line. The protein decay rate was set to (5 = 8 • lO^^.s^^. 
(B) Mean residence time ts versus mean protein level N derived from stochastic simulation (symbols) and 
our analytical approximation (lines) for four different parameter settings. Note that the analytical approxi- 
mations as well as the simulation results of the first and fourth parameter set coincide. The exponent in the 
relation cx (Na)'^ is v = 1, in accordance with equation 14. 
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1 Steady state solution of the deterministic two-stage 
toggle switch 

Using symmetric parameters a = aA = as, /? = /3a = /3b, 7 = 7A = 7b, ^ = <5a = fe, t"'" = = Tg and 
T~ = = Tg results in two following steady state solutions of equations (7), (8) and (9) in the main text: 

m^W=m^W = -^il-n) (SI) 



d^^') = dn^'^ = (S3) 



nj,W=nB^'^=^^a-v) (S2) 



1 + r? 



mj,^^)=mB^^)=-^(l + rj) (S4) 



2(3t 



nA(2)=nB(^)=-^(l + ry) (S5) 
rf^C^) = dB^^) = (S6) 



with V = y '^"st- + ^- '^^'^ ^'^^^ solution is positive, the second is negative (given all parameters are positive). 
Only the positive solution is of interest in biological systems. The monostability of this system is due to 
the monomeric protein binding. As soon as dimeric binding is introduced, the system becomes multistable, 
similar to the existing one-stage deterministic models (1-3). 

Note that for small r"*" equations (SI) and (S2) reduce to the steady state solution of a simple two stage 
expression model 

nA^^) = nB« = 

(1) (1) ^ 
7 

because 

■JOT 

for small r"*" through the Taylor approximation 

(1 + x)" « 1 + nx for |x| < 1 . 

This is expected since setting t+ to removes the interaction between both players, which will then evolve 
independently according to a two-stage expression model. 

For decreasing t~ the solution of the system approaches the origin 

UA^'^ = nn^'^ = 
rriA^^^ = TTiB^^' = 
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because rj ~ and therefore the right hand sides of equations (SI) and (S2) reduce to const • For 
decreasing r~ this term will approach 0. Hence, if proteins never unbind the promoter, the system will be 
locked forever yielding protein and mRNA levels in steady state. 

We now assess the stability of the positive solution (SI) - (S3) using standard linear stability analysis. To 
reduce the complexity of our system for the stability analysis, we apply a quasi steady state approximation 
to the DNA binding/dissociation process (c^a = = 0), reducing the dimensionality of our system to four 
equations: 

at 

— riA = /3mA - + r" (1 - il){np,)) - r+V("A)riA 
^riB = pm-Q - 5n-B + t~(1 - V(^^B)) - r+V("B)"'B , 
with il}{x) = _ ^ . The reduced system has the positive steady state solution 



UA^''^ =nn^''^ =-^{l-r,) (S8) 

with 77 = '^^g^- + 1- Notice that this is the same as the solution for mRNA and protein of the full system 
(Eq. (SI) and Eq. (S2)). We calculate the Jacobian matrix of the reduced system as 



J 



(-1 ^^•V'(nB) \ 

-7 ^-^'(nA) 

P -^ + T+V^(nA)[l + r+nA-^-^(nA)] 

V /3 -(5 + r+V'2(nB)[l + T+nB -1/^-H^b)]/ 



We evaluate the Jacobian at the steady state solution of the reduced system (Eq. (S7)-(S8)) and use the 
Hurwitz criterion to verify that all its eigenvalues have negative real part. We conclude that the system has 
one stable positive fixed point but we cannot analytically exclude the existence of limit cycles. However, 
inspection of the system's phase portrait (see Fig. S4) indicates that no limit cycles exist. 



2 Master equation 

The stochastic model of the two-stage switch is completely defined by the master equation. We define 
Pij{MA, Mb, Na, NB,t) as the probability at time t to have Ma copies of mRNAA, Mb copies of mRNAB, 
A''a copies of ProteiuA, A^b copies of ProteiuB, and the corresponding promoter configuration ij where 0(1) 
means an unbound (bound) promoter. This results in the following master equation that is split up into four 
coupled equations, corresponding to the four promoter states: 
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jVoq{Ma, Mb,Na, Nb, t) = T^E]^^roi{MA, Mb,Na, Nb, t) + t^E^^Viq{Ma, Mb,Na, NB,t) 

+ l-r^NB - t+Na + c^AiEli^ - 1) + aB{E-^^ - 1) + 1a{E+,, ~1)-Ma+ ME+j^ - 1) • Mb 
+ Pa{E^^ - 1) • Ma + /3B(i^^3 - 1) • Mb + 5a{E+^ - 1) • A^a + 5b{E+^ - 1) • Nb] 

■roo{MA,MB,NA,NB,t) 

jVll{MA,MB,NA,NB,t) = T+E+^NAVlo{MA,MB,NA,NB,t) + T+E+^NBroi{MA,MB,NA,NB,t)] 

+ [-^A - + ^a{E+^ - 1) • Ma + MEti^ - 1) • Mb 

+ f3A{E^^ - 1) • Ma + ME],^ - 1) • Mb + dA{E+^ -1)-Na + 6b{E+^ - 1) • A^b] 
•Pn(MA,MB,^A,A^B,0 

-Pio(Ma, Mb, A^a, NB,t) = T^Ej^^VniMA, Mb, A^a, NB,t) + t+ E+^NbVoo{Ma, Mb, Na, Nb, t) 

+ [-^B - ^^^A + aB(i?MB - 1) + 7a(S^, - 1) • Ma + 1b{E+^ - 1) • Mb 

+ I3a{E^^ - 1) • Ma + Pb{E],^ - 1) • Mb + 5a{E+^ - 1) • A^a + 5b{E+^ - 1) • A^b] 
•Pio(MA,MB,ArA,ArB,t) 

-Poi(Ma, Mb, A^a, A^b, i) = r^i;^^ A^AnoCMA, Mb, A^a, A^b, t) + t^E^^Vu{Ma, Mb, Na, A^b, i) 

+ [-t+ATb - + aA{E^, - 1) + 1a{E+, - 1) • Ma + TbI^^^^^ - 1) • Mb 
+ MEn^ - 1) • Ma + /3B(i^^3 - 1) • Mb + 5a{E+^ - 1) • A^a + 5b{E+^ - 1) • ATb] 

■Poi{MA,MB,NA,NB,t) . 

The shift operators E^ and increase or decrease the function argument x by one, E^ f{x) = f{x ± 1). 
To our knowledge no results have yet been published on the solution of stochastic two-stage switches. 



Reaction 
Transcription a 
Translation (3 
mRNA degradation 7 
Protein degradation S 
DNA binding t+ 
DNA dissociation t~ 



Parameter value 
0.05s-^ 

0.05 s-^mRNA-^ 
0.005 

5 • 10"-^ to 5 ■ lO-^s"^ 
1 Protein"^ 
0.1 



Table SI: Parameters of the switch model used throughout this work. Protein degradation is chosen according 
to the desired protein level n. If not mentioned otherwise, all simulations and plots are based on this set of 
parameters. 
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3 Rates 



First we derive upper boundaries for the transcription and translation rates. Transcription of DNA into 
mRNA is accomplished by the RNA-polymerase. One polymerase can process about 10-20 nucleotides (nt) 
per second in eucaryotes (4-6). As described by Alberts et al. (4) the newly elongated RNA fragment is 
immediately released from the DNA, which enables other polymerases to follow up even before the first 
mRNA has been completed. The distance d between polymerases is estimated to be around 100 nt (7). The 
rate of transcription is independent of the sequence length /, since the longer the gene, the more polymerases 
can process it in parallel. Altogether we find the maximal transcription rate a by dividing the speed of 
transcription v with the sequence length I, multiplied with the number of transcribing polymerases, 

V n 10 nt/s I „ , _i 

= T ■ ^ ~ • — = 0.1s 

I d I 100 nt 

required that enough polymerases and nucleotides are present. 

The maximal translation rate can be inferred in a similar way: Ribosomes, large complexes of proteins 
and rRNAs that translate mRNA into polypeptides, proceed with a speed v of 2 codons (= 6 nt = 2 amino 
acids) per second in eucaryotes (4). One mRNA can be processed by many ribosomes (polyribosomes) at the 
same time (4). The average space between two ribosomes is 80 nt or fa 27 amino acids (AA). (4). Therefore 
the overall translation rate for an mRNA of length n is 

I 27AA 

again independent of the mRNA length I. This corresponds to the maximally possible translation rate. The 

actual rate will be smaller when not enough ribosomes or other involved molecules (tRNA, amino acids) are 
present. We estimate the minimal translation time as l/0.074s~^ = 13.5s, which is in good agreement with 
literature, where the time needed for one translation is said to be between 20 seconds and several minutes 
(4). Notably, these transcription and translation rates only provide rough estimates of the relevant timescale. 
Throughout the manuscript, we use a transcription and translation rate oi a = j5 = 0.05 s~^, corresponding 
to an average time of 20 seconds per product, which seems to be reasonable in the context of the above 
considerations. The fact that both rates are equal is not expected to have influences on the results. 

Interactions between proteins and DNA are mediated by specific regions of the proteins, called DNA- 
binding domains, which on the one hand can recognize specific DNA sequences and on the other hand 
maintain the interaction between DNA and protein. Zinc Fingers, Leucine Zippers or Helix- Turn-Helix 
motifs are prominent examples of DNA binding domains (4). The binding between DNA and protein is 
maintained by hydrogen bonds, ionic bonds, and hydrophobic interactions. Single interactions are weak, 
but as many bonds arc formed, the binding between DNA and protein becomes stronger. The binding rates 
are very fast compared to transcription and translation processes and according to Alon (8) in the range 
of ls~^Protein~^ in Eschericia coli. The unbinding rate depends on the strength of the interaction and is 
assumed to be 10 times smaller (O.ls^^) in our model, leading to strong binding of the protein to the DNA. 
All reaction rates of the models used in this work are summarized in Table SI. 

As we showed above, the transcription and translation rates have upper bounds. The only way for a 
cellular system to further increase the abundance of proteins is to modulate the degradation rates of mRNA 
or proteins, giving longer lifetimes to mRNA and proteins. Thus, during this work we manipulate the 
degradations rates to adjust the system's protein level to a desired steady state. Notably, following the 
report of Warren et al. (9), mRNA levels are set to 10 by adjusting the decay rate. 
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4 Attractor boundaries 

In this section we assume symmetric parameters for both players (a = = ccb) • • •) for simphcity. 



4.1 Sa and 

We approximate the protein number distribution in the attractors Sa and 5b using results from Thattai and 
van Oudenaarden (10), who showed that for a simple two-state expression model, the mean and variance of 
protein numbers obey 

N = — - and a 



respectively. Thereby, we assume that in the dominating attractors, the presence of the antagonist can be 
neglected due to its marginal transcription. We define the boundary Xx of attractor using a normal 
approximation of the dominating protein's distribution as 

Xx = - Zg ■ , (S9) 

where Zg is the q% quantile of the standard normal distribution of the protein number with mean {n)x and 
standard deviation a.,,. Using q = 0.1 throughout our study assures that 99.9% percent probability mass 
of the distribution lies beyond the lower boundary. Therefore we are certain to capture all relevant protein 
numbers belonging to Sa and ^b. Using these boundaries, we can define the attractors Sa, Sb, and 5*0 
accordingly: 



Sa = {s€ S\Na >XaANb< Xb} 
Sb = {sG S\Na < Xa a JVb > Xb} 
So = {sG S\Na < Xa a JVb < Xb} 



4.2 SX and 

In the S]^/Sq states, one protein of the 'repressed' species is present, partially blocking the promoter of the 
dominating player. Therefore the dominating player is not synthesized with full rate and its protein level 
will be smaller: 

The promoter of the dominating player A is free only 100 • ^+_^^- % of the time. Therefore the average 

number of mRNAA in this attractor is m = ^+_^_^- ■ ^ (compared to ^ at full synthesis). Consequently the 
average number of ProteinA in the state 5^ will be 

iVA = m.f = ^l^."f . (SIO) 
Based on this finding one can construct boundaries in a similar way to the Sa/Sb attractors. 



5 Implications for Differentiation 

In the following, we assume that only if one player dominates the other one for a sufficiently long time td, 
downstream genes necessary for the decision process are activated (e.g. leading to chromatin remodeling), 
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and the cell differentiates. We are thus interested in the probability that the residence time T in a dominating 
attractor is longer than f^: 

LO = P{T> td) . 

This can easily be evaluated using the previous results on the residence time. Each time the switch flips, 
the system has the chance to settle in a sufficiently long residence time. What is the probability that 
after switching m times the residence time is greater than t,;? We therefore define a random variable M 
representing the number of switching events until the residence time is greater than and find that it follows 
a geometric distribution with parameter lo: 

P{M = m) = uj{l - u;)"'-^ . 

This distribution tells us how likely it is to find a cell differentiated after m switching events. It is non-trivial 
to directly relate m to the time until differentiation (further denoted by the random variable Tp) since this 
time is composed of the actual realizations of the residence time and the time the system spends in the 
intermediate attractors. As no analytic expression of the time in the intermediate attractors 5^ and Sq is 
available, we do not take into account the time spent in and 5g but only the residence time of the system 
in Sa and ^b- 

Now the time until differentiation only depends on the m — 1 independent realizations of the residence 
time (see Fig. SI A): 

m— 1 

i 

The system flips the attractors m — 1 times and resides here for a time given by until it stays in one 
dominating attractor for a time longer than td- Tp is a sum of random variables Ti but also the number of 
summands varies according to the random variable M . 

Another obstacle is the fact that the Ti are not exact geometric random variables, but truncated above 
td- For simplicity we do not make this distinction here. This is justified for small w, because uj determines 
the fraction of probability mass that is truncated. 

It can be shown that the sum Z oin geometric random variables Tj with common parameter p follows a 
negative binomial distribution: 

n 

Z = Y,Ti^NB{n,p) 

i=l 

P(Z = z)=(^^ + ^-^)(l-p)V. 

Since in our case n is drawn from a distribution itself, we have to consider the compound distribution (11) of 
M and Z to ultimately obtain the probability density function for Tp, which is a combination of how many 
switching events take place (represented by M) and the actual time the system spends in the dominating 
attractors (represented by the Ti): 
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■oo 



P{Tp = t)=Y^ P{M = m)P{Tp = t\M = m) 



m—1 




m=l 



m- 



-y 




l+p(-l+w) 



(-1 + a;)a; 



t > 



i = 



(A special case has to be considered for m = 1, where the system immediately starts differentiating. Here 
one has to define P{Tp = 0\M = 1) = oj and P{Tp > 0\M = 1) = 0, since the negative binomial distribution 
NB{n,p) is meaningless for n = 0). 

Fig. SI B shows the distribution of Tp for a mean residence time ts = Ih and a threshold of td = 5h, 
meaning that the system has to reside for more than 5h in the same attractor to promote differentiation. 
Immediate differentiation is most probable {t = 0) and the probability for higher differentiation times 
decreases quickly. Therefore the majority of cells will differentiate rapidly and only few cells will need 
longer time to differentiate. This is somewhat counterintuitive, since one would expect that it takes a long 
time until a sufficiently large residence time > td has been accomplished to promote differentiation (since 
longer residence times are less likely than shorter ones). However, due to the "odd" nature of the geometric 
distribution - where the probability for success is always the highest on the first try - the system has the 
highest chance to reach a sufficiently long residence time on the first try (refiected by the random variable 
M). Even if the system has to change attractors more often, the differentiation time will be short since the 
residence time in each attractor is most likely very small (since it is itself a geometric random variable), 
contributing only little to the overall differentiation time. 
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6 Supplementary Figures 
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Figure SI: A: Scheme of the differentiation mechanism. The system changes the dominating attractor twice 

before it settles into Sa for a time greater than t^, which induces differentiation. The residence times in 
and Sb follow a (truncated) geometric distribution. The overall time until differentiation is therefore 
(neglecting the intermediate attractors) the sum of these geometric random variables. B: Distribution of the 
differentiation time Tp for a residence time of Ih and = 5h. Small differentiation times are more probable 
than higher bigger differentiation times. Note the strong peak at < = indicating that it is most probable 
for a cell to differentiate immediately . 
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Figure S2: Quasi-potential of a toggle switch based on a one-stage expression model projected onto the 
jytotai ^totai dimensions. Notice that both axis are on logarithmic scale and are shifted by 1 in order to 
include TV^"'**' = and N^^^^ = 0. Therefore the lowest row in the plot corresponds to the case JSf^^'^^ = 0. 
The quasi-potential U{x) = — log7^'''*^(a;) is color coded where red areas reflect minima of the landscape. 
Contrary to the two-stage toggle switch model only three attractors are present. The vectors of the outflux at 
each point in state space are drawn as lines (circles correspond the the origin of the vector). Vectors are 
normalized and therefore show only the direction, not the magnitude of the field. Parameters were set in 
accordance to the two-stage switch: Protein synthesis a = 0.5 s~^, Protein degradation 7 = 5- 10~^ s~^, 
DNA binding r"*" = 1 Protein" "'^s""'^, DNA dissociation = 0.1 s^^ . 
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Figure S3: Quasi-potential of a two-stage toggle switch projected onto the and dimensions 

using a different parameter set: Motivated by Hensold et al. (12) we set the mRNA degradation rate to 
7 = 3.6 • 10~^s~^ and adjust the transcription rate a = 3.6 • 10~^s~^ to obtain the desired 10 mRNAs per 
cell (9). Furthermore we set the protein degradation rate 5 — 3.6 • 10~^s~^ as proteins are though to be one 
order of magnitude more stable than mRNAs (13). The translation rate was set to /3 = 0.0036 mRNA~^s~^ 
in accordance to (13). We observe the emergence of the same four attractors described in the main text and 
Fig. 3 therein. 
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Figure S4: Phase portrait of the deterministic two-stage toggle switch model projected onto the Na and 
TVb dimensions. The same parameters were used as in the probabilistic model of Fig. 2 and Fig. 3. The 
ODE was solved for different initial protein concentrations nj^, tlq and the mRNA concentrations where set 
to ruA = 0.01 • UA and me = 0.01 • ne, resembling the fact that for the given parameters each mRNA will 
correspond to approximately 100 proteins. Trajectories are arbitrarily colored to improve readability. Note 
that trajectories can intersect as this in only a projection of the full system. Inspection of the phase portrait 
reveals the single steady state (red dot) as predicted by Eq. (S7) - (S8). No limit cycles are visible in the 
phase portrait. 
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Figure S5: Residence time of the system in the intermediate attractors (5*^ or S^). (A) The histogram 
of the residence times in S]^ obtained from stochastic simulation suggests that the residence time in the 
intermediate attractors follows a geometric distribution. (B) Similar to Fig. 4 one observes a linear scaling 
of the mean residence time and the protein abundance. 
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